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regimes are shown to be in good agreement with established renormalization group results. We aim 
at a broad and comprehensive overview of the capabilities and shortcomings of the methods. We 
address the question to what extent the elusive low energy properties of the model are reproducible 
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I. INTRODUCTION 



In recent years both the applications for strongly 
correlated quantum impurity models and the number of 
successful approaches to harvest their physical results 
have grown enormously. Those models were introduced 
to describe the effects of magnetic transition metal 
impurities immersed in metallic hostsii^. Originally they 
were derived to capture remarkable physical properties 
like the resistance minimuni^i^ at a specific temperature 
scale Tk>^ or the anomalous magnetic susceptibility and 
specific heat of such materials. Today a whole realm of 
applications for quantum impurity models has opened. 
They describe the physics of quantum dots and wires^"— 
as well as molecular electronics^. Applications range 
from nanoelectronics all the way to quantum information 
processingiS. Their properties are essential for todays 
technological applications in single electron transistorsii 
exhibiting the Coulomb blockade effeclji^ or in devices 
dominated by RKKY interactioni^"— . The behavior of 
various magnetic phenomena and the fascinating branch 
of heavy fermion physics is described by strongly cor- 
related quantum impurity model o^^'^^ . Recent studies 
have shown that the remarkable material Graphene 
exhibits Kondo physicsi^ which may be investigated 
theoretically by virtue of quantum impurity models. 
These models have further been applied to understand 
the adsorption of atoms onto surfacesi^i"— . In addition, 
they are of theoretical importance as solvable models 
of quantum field theorie a^^i^^ . A renewed interest in 
understanding and calculating dynamic quantities of 
these models was created with the advent of dynamical 
mean- field theory (DMFT)^-^. In the foundations of 
this theory quantum impurity models have to be solved 
as an auxiliary problem. 

A wide range of methods and approximations have been 
suggested for the solution of quantum impurity models. 
They however prove to be a very delicate subject 
because standard perturbative approaches diverg 
Prominent methods to gain physical conclusions include 



a self consistent perturbative expansion^! and Bethe 
Ansatz techniques'^ for one dimensional problems. The 
low energy physics is very well described by numerical 
renormalization group (NRG)^ and in some limits also 
by functional renormalization group (FRG)^'^, and 
density matrix renormalization group (DMRG)^"— . 
There is a range of slave particle methods^i^^ available 
as well as methods based on Hubbard's X-operator 
technique^li^ and calculations using variational wave 
functions^!. Valuable physical insight has been gained 
by using equation of motion techniques applying dif- 
ferent approximation schemes^. For moderate system 
sizes the Hirsch-Fye Quantum Monte Carlo (QMC)ii 
algorithm has proven to achieve good results. In the past 
years different approaches to continuous time QMG^ 
have been applied very successfully to solve quantum 
impurity models especially in application with DMFT. 
In this context exact diagonalization (ED) methods have 
been explored to solve small systems^. 
As of today some limits of quantum impurity models are 
understood with great precision but there appear several 
gaps to be bridged. The low energy properties of these 
models are reproduced very well by renormalization 
group based approaches (i.e. NRG). These approaches 
in general have trouble to capture the high energy parts 
of the spectrum. The same may be said about QMC 
methods which if applicable yield dynamic quantities 
in imaginary time. The analytic continuation to the 
real energy axis is ill conditioned. Spectra obtained 
by, for example, the maximum-entropy method^i^ 
have a large uncertainty for higher energies. Exact 
diagonalization methods, in principle, grant access to 
low as well as high energy parts of the spectrum at the 
same time. Due to the prohibitively large Hilbert space 
however only small systems (about ten to twenty sites) 
may be treated with this method, whose low-energy 
behavior is expected to deviate from the one of the 
infinite lattice significantly. Nevertheless, the advantage 
consists in the fact that the spectral properties may be 
determined directly on the real energy axis. Besides 
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the issue of the low energy scale, also the flexibility to 
adapt to various impurity configurations and geometries 
is limited in many methods. NRG has been successfully 
applied only to the one and two impurity case so far. 
QMC approaches may suffer from the sign problem 
for more complex multiband models^. The region of 
large interaction strength is naturally difficult to treat 
in standard perturbativc/diagrammatic approaches (i.e. 
diagrammatic perturbation theory or FRG). 
In the present work we test cluster perturbation the- 
ory (CPT)^i^ and the variational cluster approach 
(VCA)i2i"— on the single impurity Anderson model^. 
The great flexibility and versatility of OPT/ VGA allows 
for obtaining approximate single-particle dynamic 
quantities and static expectation values in all parameter 
regions of any lattice impurity model with local inter- 
actions. However these many body cluster methods can 
not be expected to describe the low energy excitations 
as accurately as specifically tailored methods do. It 
is however interesting to see whether the correct low 
energy behavior may be reproduced at least to some 
extent. GPT as well as VGA bare several advantages^: 
i) They yield spectra directly on the real axis and ii) 
also the high energy incoherent part of the dynamics 
becomes available, iii) They are applicable in all pa- 
rameter regions and also at high interaction strengths, 
iv) They have the advantage of comparatively low 
computational cost for a required resolution. Our main 
goal in studying the well understood single impurity 
is to benchmark GPT/ VGA for future application to 
the not so well understood case of multiband impurity 
models in various spatial geometries. This publication 
also sets the foundations for a future extension to non 
equilibrium problems. 

The text is organized as follows. The single impurity 
Anderson model is introduced in sec. im A short review 
on GPT and VGA in this context is given in sec. IIIII 
A self consistent formulation of VGA previously intro- 
duced in the context of non equilibrium problems^ is 
presented in sec. IIII 11 Some remarks about the choice 
of variational parameters are provided in sec. IIII 21 In 
sec. IIVI we discuss the grand potential ft for infinite 
fermionic systems in relation with the VGA. Results for 
the single-particle dynamics of the SIAM are provided 
in sec. |Vl In this section also the quality of the low 
energy Kondo physics is compared to benchmarking 
results from NRG, DMRG, GT-QMG, Hartree-Fock and 
Bethe Ansatz calculations. Finally, we summarize and 
conclude our findings in sec. I VII 



II. THE SINGLE IMPURITY ANDERSON 
MODEL 

We consider the single impurity Anderson model 
(SIAM)i in real space, in one dimension 

^SIAM = ^conduction + ^impurity + ^hybridization ■ (1) 

A tight binding band of non interacting s-electrons with 
nearest neighbor {i, j) hopping is described by 

^conduction ~ "^L '^ia ~ ^ '^iu '^ja ' (2) 

1=1 cr (ij) <7 

where is the on-site energy of the particles, t is the 
overlap integral between nearest neighbor orbitals and 
i, j G {l,...,A^s} where Ns is eventually taken to be 
infinity. The operators cj^ and c,j^, respectively, create 
and annihilate electrons in orbital i with spin a. The 
impurity Hamiltonian consists of a single f-orbital with 
local Goulomb repulsion U, 

^impurity = <^f'^fl fa + U fl!^ h[ , (3) 

with creating an electron with spin a and on-site en- 
ergy e/ located at the impurity. The particle number 
operator is defined as fi^ = fl f„. Finally the coupling 
between a non interacting s-orbital and the impurity f- 
orbital is given by 

^hybridization = -V ^ c\<y f a + fl C^a > (4) 

(7 

where V is the hybridization matrix element between the 
s- and the f-orbital of the impurity atom (see fig.H] for 
illustration). 

We have set the chemical potential ^ to the center (e^) 
of the conduction electron density of states and choose 
fj, ~ Es = 0. The resonance width A is defined as 

A^nV'ps{0) = ^. (5) 

For the model defined in eq. © the local density of 
states of the conduction electrons Ps{0) is given by 
Ps(0) = In the forthcoming discussion we refer to 
the particle-hole symmetric case when we furthermore 
set ej = — All calculations are performed with t = 1 
and V = 0.3162 which yields A = 0.1. All reported 
results, except for the GT-QMG data in sec. IV Fl are for 
zero temperature. 



III. CLUSTER PERTURBATION THEORY / 
VARIATIONAL CLUSTER APPROACH 



A handle on dynamic single-particle correlations and 
expectation values is given by the single-particle Green's 
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FIG. 1: (Color online) Illustration of the single impurity An- 
derson model. The model consists of a semi-infinite chain 
of non interacting s-orbitals with nearest neighbor hopping t 
and on-site energy ts. The impurity f-orbital is subjected to 
a local on-site energy e / and local Coulomb interaction U and 
is hybridized with one of the s-orbitals (here the one at the 
beginning of the chain) via a hybridization matrix element 
V . This maps the impurity f-orbital onto site and the im- 
purity s-orbital onto site 1 in this geometry, the rest of the 
conduction band s-electron orbitals are mapped onto sites 2 
to oo. The semi-infinite non interacting chain is truncated 
at some site L. This decomposes the model into two clusters: 
an interacting cluster of variable size including the interacting 
impurity f-orbital and a semi-infinite chain of non interacting 
s-orbitals. In CPT/VCA these decomposed systems are cou- 
pled via a hopping element t. 



function Glj (w) which wc calculate within cluster per- 
turbation theory (CPT) ^^'^^ as well as the variational 
cluster approach (VCA)ii. CPT and VGA have been pre- 
viously applied inter alia to the fermionic Hubbard model 
and VGA also with great success to bosonic systems^"— 
with broken symmetry phases. The groundwork of VGA 
lies within cluster perturbation theory which is a cluster 
extension of strong coupling perturbation theory, valid to 
first order in the inter-cluster hopping. The main result 
of GPT is that the Green's function of the physical sys- 
tem G (which we call full Green's function throughout 
this text) may be obtained by a Dyson-like equation in 
matrix form 



G" 



T 



(6) 



Here g denotes the Green's function of a cluster which 
comes about by tiling the lattice of the physical sys- 
tem into smaller, numerically exactly solvable patches. 
This tiling is done by removing the hoppings between 
sites connecting such clusters. Therefore the matrix 
T = g'o"^ ^Gq"^ contains all single-particle terms connect- 
ing clusters (i.e. the inter cluster hopping which will be 
referred to as Tjnter below). The subscript o denotes the 
non interacting Green's function. To apply this approach 
to the SIAM we start by splitting the physical model un- 
der consideration (eq. (H])) into appropriate pieces. Here 
we consider a cluster decomposition consisting of two 
parts. One part, consisting of a cluster of size L, which 
contains the interacting impurity f-orbital 
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and a second, infinitely large part, the environment, 
which contains the rest of the conduction band 
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(8) 



The original Hamiltonian eq. ([TJ , defined on the semi- 
infinite lattice, may now be rewritten as 



HsiAM — ^interacting + T~L environment H~ Tj^tcr 



(9) 



Here Tintcr is the part of T describing the hopping from 
the interacting cluster to environment "cluster" , which 
is the only term not included in the two clusters. For 
the SIAM the two bare Green's functions gjntcracting ^^^^ 
gcnv (, which correspond to T^intcracting (cq. (0)) and 
T^cnvironmont (eq. ([5])), ) needed for eq. ^ may be eval- 
uated separately. This is a bit different from the usual 
application of CPT to translationally invariant systems 
which normally leads to a single cluster having discrete 
spectra, embedded in a superlattice. Therefore the ap- 
plication of CPT to this problem does obviously not suf- 
fer from issues arising due to periodization prescriptions 
for the Green's function or self-energj*^. We are deal- 
ing with two fundamentally different clusters, where one 
has a discrete (interacting cluster eq. (O) and the other 
a continuous spectrum (environment eq. (O). Due to the 
continuous spectrum of the environment a numerically 
favorable representation of the Green's function of the 
physical system G in terms of the Lehmann representa- 
tion (see for example ref. [t^ is not possible. For evaluat- 
ing quantities from the Green's function G one therefore 
has to use a direct numerical integration, which works 
best on the Matsubara axis. 

The cluster Green's function gjntoracting is determined by 
exact diagonalization of eq. ([7]). We apply the Lanczos 
algorithm^ to find the ground state and a Band Lanc- 
zos method to obtain the Green's function. The Band 
Lanczos method is initialized with a set of all annihila- 
tion and creation operators under consideration applied 
to the ground state. Thereby we obtain the so-called Q- 
matrices^ which are used to calculate the Green's func- 
tion 

ginteracting,ii(^) = ^ ^ , . \ ^31 

o V 7 / „ 

n'^t _ / 73 < 7 1 1*0 > particle part 



7=<*o|cri7> hole part 



^ I — particle part 
^ojo — uj-y hole part 

Essentially this is the Lehmann representation for zero 
temperature Green's functions. The sum over a denotes 
a sum over a possibly d-fold degenerate set of ground 
states. The sum over 7 is over a set of orthonormal 
basis-states having one particle more than the ground 
state (particle part) and one particle less than the ground 
state (hole part). 
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The Green's function of the environment g^j^^ is given an- 
alytically by the Green's function of a semi-infinite tight 
binding chain^ 

go„v,ij(^) = t^o,i-i(w) - (10) 
— i sign (3m(a;)) / oj — 



2\t\ 



+ i sign (3m(aj)) -i / 1 



2t 



where Vij is the retarded / advanced Green's function of 
the infinite tight binding chain if the infinitesimal imag- 
inary part (0"'') of uj is positive / negative. 
VGA, the variational extension of GPT, is based on 
the self-energy functional approach (SFA )^^i^^ . In the 
SFA one considers the Lcgendre transformed Luttingcr- 
Ward^ functional F[Y,], which is a universal functional 
of the self-energy, i.e. it does not depend on Go- F gen- 
erates the Green's function, i.e. 



(11) 



where /? denotes the inverse temperature. Introducing 
the (non-universal) self-energy functional 

Go] = Fm - Trln ((-G^i + S) Goo) , (12) 

(see ref. [H for a definition of Goo), one recovers Dyson's 
equation at its stationary point 



(^»[£,Go] 



-G[S] + (Go 



0. 



(13) 



Eq. (|13p is an equation for the physical self-energy S 
given the Luttinger-Ward functional Fl'S] and the free 
Green's function Go- The trace Tr is short for Tr = 
^ tr, the sum is over fermionic Matsubara frequencies 

and the small form trace tr denotes a sum over lattice 
sites and spin. The idea is that, due to its universality, 
F[S] (and thus fi[E, Gq]) can be evaluated exactly by ex- 
ploiting a different system (so called "reference system" ) 
which differs from the physical system by single-particle 
terms only. This reference system T-l' is simpler, and thus 
exactly solvable. It is defined on a cluster tiling of the 
original lattice and has the same interaction as the orig- 
inal system H. The VGA reference system is chosen to 
be a cluster decomposition of the original lattice, as the 
one introduced for GPT above. Comparing eq. for 
the full and the reference system yields 

r![E, Go] = r!'[E,go] + Tr In (- (g^ ^ - S)) 

-Trln(-(Goi-E)) , (14) 

where lowercase g denote Green's functions of the refer- 
ence system. Thus the SFA/ VGA approximation consists 
in solving eq. (|13p in a restricted range of self-energies 



E, i.e. the ones produced by the reference system. In 
this way, the space of allowed E is spanned by the set 
of single-particle parameters of the reference system, x'. 
This means that the functional ri[E,Go] feg. be- 
comes a function of those parameters 

n{x') = n'{x') + Tr In (-G(x')) - Tr In (-g(x')) , (15) 

The stationarity condition determining the physical pa- 



rameters eq. p3p is then given by 



Vx'17(x' 



(16) 



The Green's function of the physical system is obtained 
by the GPT relation eq. ©. The matrix T = gg ^ - Gq ^ 
(eq. ^) in VGA contains all single-particle terms not 
included in the reference system (i.e. Tintoi ) as well as, in 
addition, the deviation introduced by VGA, Ax = x' — x 
of the single-particle parameters of the reference system 
x' with respect to the ones of the original system x. 
In the following we fully adopt the zero temperature 
formalism, in which according expressions for the grand 
potential and related quantities may be readily evaluated. 



1. Alternative: Self consistent VGA 

In ref. [s^ we explored an alternative version of VGA 
whereby the variational parameters x' were determined 
by a suitable self consistent criterion, instead of look- 
ing for the stationary point of the grand potential 
eq. p6p . This alternative approach was introduced to 
treat systems out of equilibrium, although it can equally 
be adopted in equilibrium. The advantage of this ap- 
proach is that the solution of a self consistent equation 
is numerically easier than the search for a saddle point. 
The idea of this self consistent approach is to use a ref- 
erence system which resembles the full system best. 
The strategy is to find those values x' for the set of pa- 
rameters of the reference system which let the expecta- 
tion values of their corresponding single-particle opera- 
tors (O) cluster. x' coincide with those of the full system 
(O)cPT.x.x'- Here, the angle brackets denote expectation 
values in the cluster and the full system coupled by GPT 
or VGA respectively. Consider the on-site energies e'^ 
and e(, as variational parameters. We will look for those 
cluster parameters and e', which fulfill the relations 



cluster, e',,c', \ 'T/CPT,£/,ea,e'j,,e', 



(17) 



L-l 



L-1 



^ ("''7)cluster,e'^,£^ 5Z ("'cr ) CPT,e/ .e^ ,e'^ 



The sum is over all non interacting sites included in the 
cluster. This amounts to solving a system of non-linear 
equations in each step of the self consistency cycle. 
In general it is possible to vary each single-particle 
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parameter individually. For reasons of keeping the 
numerics tractable we use one only, which we take 
to be the same for each orbital in the chain. Extension 
to a larger number of is straightforward. To fix this 
parameter we require the average particle density on the 
non interacting sites to fulfill the condition eg . (|17l) . This 
corresponds to the condition presented in ref. |50| (see eq. 
(13) therein). In some situations, (sec below,) we will 
alternatively consider the hybridization matrix element 
V' and the intra-cluster hopping t' as variational pa- 
rameters, and proceed in an analogous way. Specifically, 
the particle number expectation values in eq. (|17l) are 
replaced by hopping expectation values. Again for t' we 
use a single variational parameter for hopping between 
all uncorrclated sites and fix it by requiring the mean 
value of hopping in the cluster and the full system to 
coincide. A discussion of this self consistency condition 
in connection with (cluster) DMFT is given in ref. [HO- 
We use an improved multidimensional Ncwton-Raphson 
algorithm to find the roots of the system eq. ([TT]). In 
some parameter regions no solution is found. 
A comparison between results obtained via the usual 
SFA, i.e. as stationary points of the grand potential 
fi, which we will now refer to as VCAq , to the ones 
obtained by the above-mentioned self consistent condi- 
tion (VCAsc ) will be given in the results section fsec.FV)). 



on-site energies we observe the grand potential ft to be 
maximal at the stationary point which is in agreement 
with results for other models. Further parameters in the 
SIAM are the hopping t and the hybridization V. As 
discussed for example in ref. l6ll . the variation of hopping 
parameters is not straightforward. For the VCAq 
approach, wc observe a maximum of fl at AV ~ —V in 
the center of two symmetric stationary points. The two 
symmetrically lying minima are equivalent due to the 
fact that the self-energy is an even power of V. As one 
tunes the parameters away from particle-hole symmetry 
this stationary point is lost in the crossover region from 
the Kondo plateau to a doubly or unoccupied impurity 
(see sec. IV Cp . In this parameter region the hopping t 
and hybridization V are probably not appropriate to be 
used as variational parameters within VCAj^. 
In the following, we always choose the set x = {V} or 
X = {V,i} for calculations at particle-holc symmetry, 
which also includes {e/, Eg}, since the variation of on-site 
energies will always yield zero deviations from the 
physical parameters and thus reproduce the CPT result 
here. For all other parameter regions it is sufficient to 
consider x = {ey, e^} as variational parameters. 



IV. GRAND POTENTIAL FOR REFERENCE 
SYSTEMS OF INFINITE SIZE 



2. Choice of variational parameters 

In VGA one can, in principle, optimize all possible 
single-particle parameters which are present in the 
original model, as well as additional ones. By adding 
bath sites not present in the original model, one in- 
cludes dynamical contributions to the cluster Green's 
function^. The numerical difficulty increases with the 
number of variational parameters. For the VGAsc case 
a multidimensional root finding algorithm has to be 
adopted. For the VGAfj case, a saddle point in many 
dimensions has to be located. Since the allowed set of 
variational parameters limits the search space for the 
self-energies one will find a solution in this restricted 
space only. It is therefore desirable to vary as many 
single-particle quantities as possible. A balance has to 
be found between a large space of available self energies 
and numerically feasible multidimensional algorithms. 
Many works have addressed the question of which 
parameters are the most important to vary and how the 
choice of variational parameters will influence or limit 
the results^. As discussed in refs. [sqUboI it is important 
to include an overall chemical potential as a variational 
parameter in order to preserve thermodynamic consis- 
tency. As a compromise, we will take two variational 
parameters x = {e/, Cs}, which cover the overall chemical 
potential. Note that this amounts to shifting an overall 
on-site energy in the whole cluster plus an extra inde- 
pendent shift at the correlated site. For the variation of 



The reference system consists of two parts, a finite in- 
teracting system and a non interacting system of infinite 
size, the environment. J7'(x') is given by the sum of the 
grand potentials of the interacting cluster (^^interacting) 
and of the environment (J^env) d which correspond to 
"Hintoracting (eq. ©) and ?ionvironmcnt (eq. ©)). Here we 
outline how to determine the grand potential for such 
kinds of reference systems. For the Green's function 
G within the GPT/VCA approximation the Dyson equa- 
tion is given in eq. The Green's function and T have 
the block structure 



G 



Gcc Gc 

Gpr Gr 



T T 

T.r. 



Up to this point all matrices involving environment in- 
dices have infinite size. As far as the Green's function 
itself is concerned this is no problem as we arc primarily 
interested in Gcc for which the Dyson equation reduces 
to 



icc^ceGec 



Gcc — gcc 
Gee — ^r^ee^GcGcc ; 

and therefore 

Gcc = gcc + Scc^ccGcc 
^cc • Tec TcegecTec ■ 

A bit more tedious is the elimination of the explicit 
dependence on the environment part of G, as far as 
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the grand potential eq. ^TE\\ is concerned. We start out 
from a form of the grand potential functional given by 
Senechal^ 



An:=n-n' = -Tr In 11 - Tg 



(18) 



In app.|A]it is shown that Ail can be expressed solely in 
terms of cluster quantities 



An 



-Tr In 1 



(19) 



Along the lines outlined in ref. |58|, the resulting integral 
can be regularized and expressed as 



O — O' — O' 

'^'cnv ^'interacting 



tr(T) 



-E 



do; In 



det a 



(20) 



The quantities ft' are the grand potentials of the uncou- 
pled reference system. The constant infinite contribution 
n'^nv is absorbed into the definition of il. It plays no 
further role as it does not depend on the variational 
parameters. This integral may be evaluated as suggested 
in ref. m by integrating from to Ai, from Ai to A2 and 
from A2 to cx). Ai and A2 represent two characteristic 
scales in the problem (for example the smallest /largest 
eigenvalue of the Hamiltonian matrix). For the last 
part of the integral a substitution w = ^ is performed. 
We use an adaptive Gauss Legendrc integrator for the 
evaluation. 



V. RESULTS 

We have evaluated several benchmarking dynamic 
quantities of the SIAM. In the following, results for the 
impurity density of states will be presented and com- 
pared to ED, NRG and DMRG data. We will elaborate 
on the strengths and weaknesses of the methods as well 
as the comparison of CPT to VGA. Furthermore, we 
will discuss the relation between VGAsc, where the 
variational parameters are determined self consistently 
via eq. (jl7p and VG An , where the variational param- 
eters are defined at the stationary point of the grand 
potential. We will show that the Kondo resonance is 
reproduced within the framework of GPT/VGA and that 
the variational results fulfill certain analytic relations 
like the Friedel sum rule (eq. ((23)) ). The method will 
be shown to provide reasonably accurate results in a 
wide range of parameter regimes of the model. Low 
energy properties related to the Kondo temperature Tk 
will be discussed in context with renormalization group 
results. The imaginary frequency Green's function and 
self-energy will be compared to GT-QMG results. 



A. Even-Odd Effect - choice of the impurity 
position 

GPT/VGA rely on the Green's function of an in- 
teracting cluster of size L which is obtained by exact 
diagonalization. Due to this fact, it is unavoidable 
that some effects of the finite size cluster influence 
the solution of the full system. (Except in the case of 
vanishing interaction strength, i.e. U = 0.) Therefore 
suitable clusters have to be chosen on a basis of physical 
results. Some aspects of this are discussed by Balzer 
et al^ in the context of DMFT and VGA and by 
Hand et al^ in the context of DMRG. In this work 
we consider interacting clusters of even size only. For 
these systems the ground state does in general not suffer 
from spin degeneracy. Furthermore, the spatial position 
of the impurity is important. This can be inferred 
from the bath's density of states, which vanishes for 
a; = at every second site. It may also be seen in the 
structure of the ground state, for which the size of the 
degenerate sectors alternates with the geometrical size 
of the cluster. Throughout this work we position the 
impurity f-orbital at the beginning of the infinite chain, 
although essentially the same results are achieved by 
attaching it to an s-orbital at site two, four, etc. inside 
the chain. 



B. Spectral properties 

The single-particle spectral function is obtained 
from the retarded Green's function G^'™* (see e.g. 

I I lb \ "-J 

ref. ^ 



1 

-— ym 

TT 



(^)) 



(21) 



The diagonal element at the impurity f-orbital ^^^(w) 
describes the impurity density of states pj(w). A 
physical property of the SIAM which poses a challenge 
to numerical methods is the Kondo-Abrikosov-Suhl 
resonance often referred to as Kondo peak^. It arises 
in the parameter regime where the magnetic moment 
of the impurity is screened by the conduction electrons 
to form a singlet stated. The particle-hole symmetric 
model lies in the center of this Kondo region. This quasi 
particle excitation is for example not captured in mean 
field approaches. With increasing interaction strength U 
the numerical solution becomes increasingly challenging. 
In this section we elaborate on the results for the density 
of states in the particle-hole symmetric case. 
A comparison of the local single-particle spectral func- 
tion at the impurity f-orbital as obtained by exact 
diagonalization and VGA^ is shown in fig.Hl The ED 
result is for a ten site system with open boundary 
conditions. The VGAq result is for an infinite reference 
system, where the interacting part of the reference 
system was taken to be of size L = 10. An ED treatment 
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FIG. 2: (Color online) Comparison of the local single-particle 
spectral function at the impurity f-orbital at particle-hole 
symmetry as obtained by exact diagonalization (ED) of a 
ten site chain (black) and VCAn (cyan). VCAn was used 
with two variational parameters: the hopping t and the hy- 
bridization V considering a length of the interacting part of 
the reference system of L = 10. All data shown is for an inter- 
action strength of U/A = 12. All results have been obtained 
for a large numerical broadening 0^ = 0.05. The inset shows 
a zoom to the low energy region. 



of a finite-size SIAM can not reproduce the low energy 
resonance at zero energy in the single-particle spectral 
function (see also app.[B|). It consists, in the particle- hole 
symmetric case, for an even number of sites (and open 
boundary conditions), of symmetrically lying excitations 
which shift closer to zero for increasing system size and 
represent a large energy scale. For an odd number of 
sites a pole in the local single-particle Green's function 
of the impurity f-orbtial will be present at w = 0. CPT 
as well as VGA are able to reproduce finite spectral 
weight at cj = (even for 0+ 0), since these methods 
are formulated for an infinite system. The finite-size 
structure in the high energy incoherent part of the 
spectrum, owing from the excitations of the interacting 
part of the reference system, is strongly reduced in 
VGAn . 

Results for the single-particle spectral function eq. (PT|) 
of the impurity f-orbital are shown in fig. [3] for four 
different interaction strengths U/A = 4,8, 12 and 20. 
As a reference, the spectra obtained with NRG and 
DMRG from Peters^ are plotted. Renormalization 
group approaches like NRG are especially suited to 
reproduce the low energy quasi particle excitations 
of this model and therefore serve as a reference for 
our data. The spectra of Peters were obtained for a 
fiat conduction electron density of states, which was 
mapped by linear discretization in energy space onto 
the corresponding orbitals of a semi-infinite chain. Our 
model is based on a semi-circular density of states 
of the conduction electrons. The low energy part of 



the spectra is comparable because we have chosen the 
only relevant parameter for the low energy part of 
the spectrum: A accordingly. This parameter fully 
determines the influence of the conduction electrons on 
the impurity f-orbital for low energies and therefore the 
low energy part of the spectrum. The high energy part 
of the spectrum may deviate slightly and is not directly 
comparable but yields a crude reference. In addition, 
we have chosen here a very large numerical broadening 
of 0+ = 0.05 for reasons of comparison only. This value 
was used in the DMRG calculations and is needed there 
to obtain spectra using the correction vector method. 
This influences the width and the height of the Kondo 
resonance, located at ui = 0. The CPT spectral weight 
at oj = appears too broad in the plot in comparison 
with the NRG result. This is only partly due to a large 
numerical broadening. Due to the nature of the GPT 
method we cannot expect it to reproduce the low energy 
spectrum as well as RG calculations do. The height of 
the Kondo resonance appears too small in this flgure 
because of the large 0+. It converges with 0"*" — ?> 10~^ 
to the result predicted by scattering theory (see fig. [5] 
and sec. IV C|) . The high energy incoherent parts of the 
spectrum located at uj ~ — e/ and lo ~ —ef + U develop 
more and more with increasing length of the interacting 
part of the reference system L. A comparison of the 
center of gravity of the high energy incoherent part of 
the spectrum of the i = 14 site GPT result and the fifty 
site DMRG result is in reasonable agreement. There are 
spurious structures in the spectral density, originating 
from the cluster Green's function of the finite system, 
preventing continuous spectra to form. We would like 
to note that the accurate determination of the Green's 
function of the reference system is of prime importance. 
An inaccuracy in pole-positions or pole-weights for very 
small but non-vanishing weights will yield spurious 
artifacts in the spectra in the vicinity of a; = 0. 
To improve on the result of CPT we considered the 
hopping matrix element t and the hybridization matrix 
element V as variational parameters. The parameters 
used for the evaluation of the reference system were 
determined with two different methods. VCAn results 
are depicted in the plot for a length of the interacting 
cluster of L = 10. As shown in the figure this method 
strongly reduces the finite size peaks in the high energy 
incoherent part of the spectrum. The width of that part 
of the spectrum is reproduced correctly for high values 
of U where the full width at half maximum (FWHM) 
within VGA is given by « 1.9 A. This comes very 
close to the expected 2 A^^i^^ of the high energy atomic 
excitations. VGA improves the spectral properties of 
the Kondo resonance with respect to GPT, bringing it 
closer to the fifty site DMRG result. The data obtained 
using the self consistent VGA approach VGAsc agree 
very well with the result based on VGAn on the position 
of the spectral features. The respective weight however 
disagrees for low values of interaction strength U, which 
is due to the different values predicted for the variational 
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parameters by the two procedures. One should note that 
the two broad Lorentzian high energy peaks (in VCAfj 
as weU as VCAsc ) consist of many excitations which 
wiU be revealed upon repeating this calculation with 
smaller 0+. 

A more detailed look on the spectral region of the 
Kondo resonance is provided in fig. El The CPT/VCAa 
data is compared to NRG and FRG data as well as 
results obtained from a restricted Hartree-Fock calcu- 
lation from Karrasch et al^. The CPT/VCA resuhs 
are plotted for lengths of the interacting part of the 
reference system L = 2,4,6,8 and 10 for two different 
sets of parameters. The results for higher L are always 
located towards the center of the figure. The results 
corresponding to the resonance at w = were obtained 
for the particle-hole symmetric model. For this set of 
parameters we used the hybridization y as a variational 
parameter. The second peak shown centered around 
w/A « 0.8 corresponds to a parameter set right at the 
border of the Kondo region. The variational parameters 
used away from particle-hole symmetry are x = {e/, Cs}- 
One can see that the OPT result is not converged for 
L = 10 site interacting clusters yet. In contrast, the 
VCAn result seems to converge much faster. Although 
in the plot it looks like the VGA result does not improve 
much upon a restricted Hartree Fock calculation, we 
will show in the following that GPT/VCA yields results 
in all parameter regimes of the SIAM which cannot be 
reproduced within a mean field treatment. 

The variational parameters obtained for the two sets 
of parameters shown in fig.[S] are presented in fig. El In 
addition to the VGAq parameters, which were used for 
the results above, the variational parameters obtained in 
VGAsc are also depicted. We plotted the difference of 
the parameter of the reference system x' to the physical 
parameter x: Ax. All deviations Ax appear to converge 
to zero with increasing length L of the interacting part 
of the reference system. Notice that the self consistent 
approach always leads to a Ax of greater magnitude with 
respect to VCAn . Remarkably, the spectra obtained by 
VGAn and VGAsc for the parameter set x = {e/,es} 
are in very good agreement even though the variational 
parameters are rather different. The most striking 
difference is that the self consistent approach yields a 
negative Acf while the based VGA yields a positive 
Ae/. This is however compensated by the different Ae^. 
Using the hybridization 1/ as a variational parameter, 
the AV obtained by VGAo and VGAsc agree rather 
well. Remarkably, the resulting density of states is very 
different, which shows that the calculation is extremely 
sensitive to this parameter. 

A spatially resolved image of the spectral function, 
calculated with GPT, for the parameter set used in 
fig-El (c) is shown in fig.lH The qualitative picture 
would be the same in VGA, merely the structures are 
slightly shifted. This view reveals how the perturbation, 
introduced by the impurity, is fading away slowly in 
an alternating fashion. At every second site away from 



the impurity a dip at w = is present, which is usually 
referred to as Fano dip. 

The low computational effort of GPT/ VGA proofs ad- 
vantageous for calculating spectra. The VGA procedure 
(for a twelve site interacting cluster) usually converges 
in minutes to hours on a standard workstation PG, while 
more demanding numerical methods often need days 
to a week to converge. Furthermore, the spectra are 
exactly determined from the Lehmann representation 
and no ill-posed analytical continuation is required in 
comparison to methods working in imaginary time or 
imaginary frequency space. To our knowledge the most 
accurate spectra available for this model so far are 
published in ref. IH. 

Overall one may conclude that GPT, YCAq and VGAsc 
reproduce a Kondo resonance, which fulfills the Friedel 
sum rule (cq. ^) for 0+ 10"^. The VGA resuhs 
improve drastically upon the GPT data which may be 
seen in a much faster convergence in L and a suppression 
of finite size effects especially in the high energy part of 
the spectrum which in addition has the expected width 
within VGA. VGAn and VGAsc agree rather well on the 
position of the spectral features. However they assign 
very different spectral weight to them at low values of 
interaction strength U. 



C. Impurity density of states and occupation 

The occupation of the impurity f-orbital is given at 

r = oby 

1 1 f°° 
+ - dw^e{G^ff{iLo)) . (22) 

This integral may be evaluated from the imaginary fre- 
quency Green's function, which in turn is directly acces- 
sible within GPT/VGA. 

To see whether GPT/VGA are good approximations in 
all parameter regions of the SIAM, we vary the on-site 
energy of the impurity e/ at fixed interaction strength 
U . The local impurity density of states at the chemical 
potential (w = = 0) and the impurity occupation num- 
ber are plotted for various lengths of the interacting part 
of the reference system L = 2, 4, 6 and 8 for the same 
model parameters. The VGAn result is shown in fig. [3 a 
VGAsc calculation in fig.Eland the GPT data in fig.El 
We start out by discussing the VGAn result (fig. [7]). The 
variational parameters x used within VGAn are the on- 
site energy of the impurity e f and the on-site energies of 
the uncorrclatcd cluster sites eg. The density of states 
p/(0) displays a pronounced plateau which is related to 
the existence of a quasi particle peak (Kondo resonance) 
pinned at the chemical potential. The parameter regions 
leading to an empty (— £/ < 0) or to a doubly occupied 
(— ey > U) impurity do not show a pinning of the Kondo 
resonance at the Fermi energy, as expected. In the half 
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FIG. 3: (Color online) Single-particle spectral function at the impurity f-orbital at particle-hole symmetry for different interac- 
tion strengths U . The interaction strengths shown are C//A = 4 in the upper left figure (a), [//A = 8 in the upper right figure 
(b), [//A = 12 in the lower right figure (c) and f7/A = 20 in the lower right figure (d). Each plot shows the results obtained by 
CPT for a length of the interacting part of the reference system of L = 14 (dashed-red), VCAn with two variational parameters: 
the hopping t and the hybridization V which are determined by the stationary point of the grand potential $1 at a length of 
the interacting part of the reference system of L = 10 (blue), VCAgc with the same variational parameters determined self 
consistently at a length of the interacting part of the reference system of L = 10 (cyan). All results have been obtained for 
a large numerical broadening 0"*" = 0.05. As a reference the NRG and DMRG results of Peters^ are plotted in yellow and 
dash-dotted-dark brown respectively. 



filled region which lies in between, virtual spin fluctu- 
ations lead to a pronounced quasi particle peak at the 
chemical potential. We observe that the result converges 
with increasing length of the interacting part of the ref- 
erence system L to the physically expected result. Due 
to the variational parameters considered, the deviations 
of the results as a function of L are rather small as com- 
pared to CPT where the results change signiflcantly with 
increasing size of the reference system (see flg-IH])- We ex- 
pect CPT calculations in the empty or doubly occupied 
regions to converge rather fast (within a few sites) while 



calculations in the Kondo regime, and particularly in the 
crossover region, may fully converge only at very large 
(i.e. exponentially) sizes of the reference system^. This 
is inferred from the spin-spin correlation function in the 
cluster which is observed to decay sufficiently fast out- 
side the Kondo plateau (i.e. it is effectively zero at the 
boundary of the cluster) but shows long range correla- 
tions inside the plateau. The VCAsc results are obtained 
with one variational parameter x = {e/}. The reason for 
not using x = {£^,6^} again is that the result is almost 
the same as the one obtained with VCAji (see fig-E]). 
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FIG. 4: (Color online) The local density of states eq. ^21^ is 
shown resolved in real space. The spectrum was obtained 
using CPT on a L = 14 site interacting cluster. The impurity 
parameters were U/A = 12, e//A = —6 and the numerical 
broadening was set to 0"*" = 0.05. The spectrum shown in fig. [3] 
(c) corresponds to the data shown for the impurity f-orbital 
located at site in the plot. The density plot is shown with 
a logarithmically scaled coloring from blue indicating zero to 
red indicating high values. 



However in some (small) parameter regions the numeri- 
cal evaluation becomes difficult. The VCAsc data shown 
in fig. |8] shows a clear improvement as compared to CPT 
but does not reach the quality of the VCAn result in 
terms of convergence in system size. 
The Friedcl sum rule (FSR)^iSi^ provides an exact re- 
lation between the extra states induced below the Fermi 
energy by a scattering center and the scattering phase 
shift. It also holds true for interacting systems. This 
gives a relation between the f-orbital occupation {n-^}, 
and the density of states at the Fermi energy: 



P/(0) 



1 
ttA 



(23) 



In our case the mean occupation in the Kondo regime 
is (n^) « 1. Since both, the occupancy of the f-orbital 
and the magnitude of the local density of states at 
the Fermi energy, can be evaluated independently, we 
can check the validity of the Friedel sum rule in our 
approximation. Results are shown in fig. [7] applied to 
the L = 8 site VCAn data. The VCA^ results fulfill the 
Friedel sum rule almost in the whole Kondo region. At 
the crossover to an empty or doubly occupied impurity 
the Friedel sum rule is not fulfilled exactly any more but 
approximated very well. Further outside the agreement 
is again excellent. The variational parameters of VCA 
are crucial to fulfill the Friedel sum rule as can be 
seen from a CPT calculation (fig. [9]) which violates it 
in all parameter regions. It appears that VCA^ with 
variational parameters x = {e/,es} naturally drives 



FIG. 5: (Color online) Magnification of the Kondo resonance 
in the density of states of the impurity f-orbital. Shown are 
calculations for two different sets of parameters. The reso- 
nance at oj = corresponds to the particle-hole symmetric 
model: U/A = 20, e//A = —10, while the resonance away 
from zero corresponds to a set of parameters right at the edge 
of the Kondo region: U/A = 20, e//A = 0. For compari- 
son we show NRG (yellow) and FRG (dark brown) data as 
well as results obtained from a restricted Hartree-Fock cal- 
culation (blue) from Karrasch et al^ (The NRG results are 
partially hidden by the FRG results.). The CPT result (cyan) 
is shown for lengths of the interacting part of the reference 
system L = 2, 4, 6, 8 and 10. Results for higher L are always 
located towards the center of the plot. In the particle-hole 
symmetric case VCAn (magenta) was performed with varia- 
tional parameters x = {V} for L — 2, 4, 6, 8 and 10. Away 
from particle-hole symmetry VCAn was performed with vari- 
ational parameters x = {e/, £3} for the same lengths of the in- 
teracting part of the reference system L. For the CPT/VCA 
calculations a numerical broadening of 0^ — 10~^ was used. 
The inset shows a zoom to the top region of the peaks. 



the system to fulfill this condition. The VCAsc result 
(fig-Hi) violates the sum rule too. This is not a feature of 
VCAgc in general but has to do rather with the choice of 
variational parameters, which was just x = {e/} in this 
case. The VCAgc result for two variational parameters 
x = {ef,es} looks qualitatively like the respective VCAq 
result. 

Scanning the interaction strength U at fixed im- 
purity on-site energy ef confirms the presence of the 
Kondo behavior. Shown in fig.[in] arc results obtained 
with VCAfj using the same variational parameters 
x = {ef,es} as above. In the weakly correlated part 
{U /A < 5) the density of states at the chemical potential 
is low. The intermediate region (5 < J7/A < 15) signals 
the crossover to the Kondo regime. For larger U the 
Kondo regime is reached with an impurity occupation 
of (n^) w 1, which may be inferred from the Friedel 
sum rule. In the inset of the figure, the CPT results for 
the same lengths of the interacting part of the reference 
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FIG. 6: (Color online) Evolution of the variational parame- 
ters for the data shown in fig.[5] Shown is the difference of the 
parameters of the reference system x' to the physical param- 
eter x: Ax = x' — X. Parameters obtained by VCAn (crosses) 
are compared to those obtained by VCAsc (circles). The 
variational parameters Ae/ (dark brown) and Ae^ (yellow) 
correspond to the calculation away from particle-hole sym- 
metry in fig.[S] while the parameter AV (olive) corresponds 
to the calculation at particle-hole symmetry. Lines are only 
guides to the eye. 



system L are shown. The CPT results are by far not 
converged for the interacting cluster sizes considered 
here. This emphasizes the importance of the variational 
parameters. 

Our results in fig. [7] and fig.fTU] agree very well with 
those of calculations based on X-operator technique 
exercised by Lobo et ali^. In their work a strong 
coupling perturbation theory is applied starting from 
the Anderson molecule as a basis and using the Friedel 
sum rule as a condition to fix the position of an infinitely 
narrow conduction band. 

Analytic considerations (see app. |B|) allow insight into 
the behavior of the Friedel sum rule in ED, CPT and 
VGA. There it is shown that ED always has to violate 
the Friedel sum rule while CPT always fulfills it in the 
particle-hole symmetric case. This comes about in the 
first place because the height of the Kondo resonance at 
w = does not depend on the self-energy. A pinning of 
the Kondo resonance however can only be achieved via 
the improved self-energy contributions obtained within 
VCA. 

The results of this section clearly show that VCA is 
able to capture the basic physics of the SIAM in every 
parameter region. The improvement obtained by going 
over from CPT to VCA is crucial to fulfill exact analytic 
relations. Moreover we have shown that CPT/VCA is 
clearly superior to ED calculations. VCA is, even at 
small L, capable of fulfilling the FSR also away from 
particle-hole symmetry due to a pinning of the Kondo 
resonance at the Fermi energy. This pinning can be 
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FIG. 7: (Color online) Density of states of the impurity f- 
orbital (solid lines) obtained via VCAn at co — and av- 
erage occupation of the impurity (dashed lines) for differ- 
ent lengths of the interacting part of the reference system 
L — 2, 4, 6 and 8 (blue, green, red and cyan) as a function 
of the impurity on-site energy e / . The Coulomb interaction 
U is kept constant at U /A — 20. The numerical broaden- 
ing used is 0^ = 10"*^. The set of single-particle parameters 
considered for variation within VCAn is x — {e/,es}. Note 
that here the point £/ = — y corresponds to the particle-hole 
symmetric case. The Friedel sum rule (eq. (|23|l 'l was applied 
to the L — 8 result (dotted- violet). It is fulfilled to a very 
good approximation in the Kondo region and far outside of 
it. Small deviations from the Friedel sum rule arise at the 
crossover region to an empty or doubly occupied impurity. 
The inset shows a zoom to the Kondo plateau. 



attributed to the better approximation of the self-energy 
of VCA with respect to CPT. 



D. Crossover diagram 

To delve into the CPT/VCA results for the whole 
parameter range of the SIAM a "phase diagram" is 
presented in this section. This should be understood to 
be a mere scan of the parameters U and e/ because the 
model does not undergo a phase transition. The density 
of states of the impurity at the chemical potential 
Pf{0) is shown in fig.[TT] in a density plot. This figure 
essentially shows the height of the Kondo resonance as 
a function of interaction strength and on-site energy of 
the impurity. The different regimes of the SIAM, as 
obtained by an atomic limit calculation, are indicated as 
black lines. These lines divide the physics into regions 
where the impurity is doubly, singly or not occupied. In 
the singly occupied region ( ^ > |e/ + ^|) local moments 
and their screening arc expected to appear. This region, 
which bestrides the cone enclosed by black lines, is the 
region where Kondo physics may take place within this 
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FIG. 8: (Color online) Density of states of the impurity f- 
orbital (solid lines) obtained via VCAsc at w = and av- 
erage occupation of the impurity (dashed lines) for differ- 
ent lengths of the interacting part of the reference system 
L = 2,4,6 and 8 (blue, green, red and cyan) as a function 
of the impurity on-site energy e/. The Coulomb interaction 
U is kept constant at (7/A = 20. The numerical broadening 
used is 0+ = 10"". The set of single-particle parameters con- 
sidered for variation within VCAsc is x = {e/}. Note that 
here the point £/ = — ^ corresponds to the particle-hole sym- 
metric case. The Friedel sum rule (eq. (|23p l was applied to 
the L = 8 result (dotted- violet). It is fulfilled in a region of 
nf < 0.4 and n/ > 1.6. 
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FIG. 9: (Color online) Density of states of the impurity f- 
orbital (solid lines) obtained via CPT at a; = and av- 
erage occupation of the impurity (dashed lines) for differ- 
ent lengths of the interacting part of the reference system 
L = 2, 4, 6 and 8 (blue, green, red and cyan) as a function 
of the impurity on-site energy e/. The Coulomb interaction 
U is kept constant at U /A — 20. The numerical broaden- 
ing used is 0^ = 10"". Note that here the point = — ^ 
corresponds to the particle- hole symmetric case. The Friedel 
sum rule (eq. (|23[l ) was applied to the L = 8 result (dotted 
violet). It is drastically violated. However the results are far 
from converged for the small lengths of the interacting part 
of the reference system considered here. 



approximation. The parameter regions where the im- 
purity is empty or doubly occupied lie above and below 
this cone. More sophisticated methods will lead to a 
smearing out of the border of these regions and introduce 
a crossover area with competing effects. A boundary 
expected between a single resonance and a spurious 
local moment behavior where the single resonance is 
split into two for spin up and spin down respectively 
is obtained by mean field thcorj*^. In the mean field 
approach the density-density interaction of the impurity 
Hamiltonian is replaced by a spin dependent density 
exposed to the mean contribution of the other spins' 
density. The mean field boundary is then obtained by 
replacing the mean field parameters (n^) and (n^) by 
the particle number n and magnetization m. Setting the 
magnetization to m = 0^ in the self consistent equations 
yields the implicit result Uc ~ 7rA(l + cot(7r-^^^)^) and 
6/,, ^ A(cot(^iM) + 1(1 _ (^nf)) (1 + cot(^i^)2)). 
The plot shows that the Kondo plateau is reproduced 
very well by VCAj^. The results appear almost con- 
verged for lengths L sa 6 of the interacting part of the 
reference system. Increasing L yields better results in 
the crossover region. Results obtained by means of CPT 
do not reproduce the Kondo plateau very well for small 
L. 

The average impurity occupation for the same pa- 
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FIG. 10: (Color online) Density of states of the impurity 
f-orbital at = for different lengths of the interacting 
part of the reference system L = 2, 4, 6 and 8 (dark brown, 
cyan, olive and magenta) as a function of the interaction 
strength U. The impurity on-site energy e/ is kept constant 
at 6/ /A — —10. The numerical broadening is chosen to be 
0+ = 10"". The set of single-particle parameters considered 
for variation within VCAn is x = {e/,es}. The inset shows 
the CPT results. 




FIG. 11: (Color online) In this plot a "phase diagram" of the 
SIAM is shown. The quantity on the z-axis is the density 
of states in the impurity pf at oj = (i.e. the height of 
the Kondo resonance). The results are obtained with VCAn 
for a set of variational parameters x = {e/,es}, L = 6 and 
0^ — 10~®. The black line indicates the different regions 
obtained from an atomic limit calculation. In the right cone 
local moments are to be expected. While in the upper region 
the impurity is expected to be empty and in the lower half 
to be doubly occupied. The blue curve shows the onset of a 
spurious magnetic state as obtained by a mean field treatment 
(see text). 



ramctcr region is shown in fig. 1121 The result obtained 
with VCAn clearly shows the Kondo plateau where the 
impurity is singly occupied. The parameter regions of 
a doubly occupied or empty impurity lead to a density 
of states in the impurity which is zero at the chemical 
potential (compare to fig.[TT|). 

The results of this section have been obtained using 
VCAn with variational parameters x — {e^,es}. It 
should be noted that using only x = {e/} already yields 
good results. As mentioned in sec. IV C| CPT needs very 
large values L to yield the same quality of the results as 
VCA does with much smaller values of L. 



E. Low energy properties, Kondo Temperature 

In this section we examine the low energy properties 
of the symmetric SIAM. In the strong coupling limit a 
single scale, the Kondo temperature Tk, governs the low 
energy physics^. 

The Kondo temperature Tk is known from Bethe Ansatz 
results for the particle-hole symmetric SIA M^^'^° 
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FIG. 12: (Color online) Average particle density in the impu- 
rity {n^) as a function of U and e/ for the same parameters 
as in fig. llll 



This scale, which is inversely proportional to the spin- 
flip rate of the impurity, divides the physics of the SIAM 
into two regions: A local moment behavior of the impu- 
rity, where the spin is free, and a low temperature region 
where the local spin and the conduction electrons become 
entangled and form a singlet stated. 
Quantities which depend inversely on Tk are the effec- 
tive mass m* and the static spin susceptibility Xm- The 
Kondo temperature may furthermore be extracted from 
the width or weight of the Kondo resonance in the local 
density of states of the impurity f-orbital. We investi- 
gate and compare the results for the scale Tk obtained 
from the direct determination of Tk (from the FWHM 
and the spectral weight of the Kondo resonance) and the 
inverse quantities m* and Xm- We find that the results 
of all four measurements turn out to yield the correct 
qualitative behavior in VCAn . However, in a region 
where the dependence of Tk is exponentially dependent 
on the interaction strength U the exponential prefactor 
is not predicted correctly. Therefore we introduce a scal- 
ing factor 7 (eq. (HI])) which turns out to be the same for 
all four ways of determining Tk- In particular this fac- 
tor is independent of the set of model parameters used. 
The scaling factor may be calculated semi-analytically 
for a reference system consisting of a two site interacting 
cluster and the semi-infinite environment within VCAn 
and VCAsc (x = {^})- The calculation for VCAn leads 
to an integral expression for the stationary point of the 
grand potential with respect to AF from which the 
optimal /SV can be obtained numerically (see app. [Cjl . 
The Kondo scale may be determined from the so obtained 
values of V'{U) = IW[U) + V hy 

TK{U)^{pp-^\ (25) 
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This leads to a perfect exponential behavior as defined 
in eq. (IMl) with 

7 = 0.6511 . 

The issue of obtaining an exponential scale but not the 
correct exponent for the functional dependence on U is 
common to various approximate methods (for example 
variational wave functions where the issue was cured 
by introducing an extended Ansatz by Schonhammer— , 
saddle-point approximations of a functional integral 
approacb^^ or FRGl^). A faint analogy may be drawn 
here to Gutzwiller approximation, where an exponential 
energy scale in U arises by a renormalized hybridization 
parameter V—, which is also the case for VCAq . 
The self consistent calculation for VCAsc also leads 
to an integral expression for the determination of AV. 
This expression is obtained by requiring the expectation 
values of the hopping from the impurity f-orbital to the 
neighboring site in the reference system to be the same 
as the expectation value in the physical system. This 
procedure does not yield an exponential scale in U. The 
optimal cluster parameter V' shows spurious behavior 
as a function of U. Wc conclude that VCAsc with 
X = {V} cannot reproduce the low energy properties 
of the SIAM even qualitatively, while VCAq yields the 
correct behavior apart from an exponential factor. 



Effective mass - quasi particle renormalization 



The effective mass m* 
renormalization^ 



is defined as the quasi particle 



m*{U) _ ^ d[SmS^^(»u;,[/)] 
TO*(0) duj 

d[3mG^^(iw,[/)] 



=0+ 



d[3?mGy^(iw,0)] 



j=0+ 



j=0+ 



(26) 



where we introduced the dependence on the interaction 
strength U explicitly. In the Kondo regime, this quantity 
becomes inversely proportional to the Kondo tempera- 
ture. 

We want to answer the question whether the Kondo scale 
is approximately captured by CPT/VCA or not. There- 
fore we compare the functional form and the exponent 
obtained from the effective mass and the analytic result 
for Tk, cc^. ((24)) . The result for the effective mass ob- 
tained within VCAo is shown in fig. 1131 The variational 
parameter used was x = {V}. The functional form is re- 
produced well by VCAo (i.e. it starts out quadratically 
and goes over to an exponential behavior in the Kondo 
region). However the exponent (g^) is not reproduced 
correctly. VCAq yields a lower exponent of (Tss)- 
The factor 7 is defined in app. [Cl determined from a 



3 





NRG 




FRG 




VGA^, L=6 




CPT, L->>= 


+ 


CPT, L=2 





CPT, L=4 


A 


CPT, L=6 




CPT, L=8 


□ 


CPT, L=10 




U/A 



FIG. 13: (Color online) Effective mass m* of the Kondo reso- 
nance eq. (|26p as a function of interaction strength U . We plot 
CPT results for lengths of the interacting part of the reference 
system L = 2, 4, 6, 8 and 10 (magenta), the to L — >■ 00 extrap- 
olated CPT result (olive) , as well as VCAn results (blue) . The 
data points for the CPT result in the low U region are not 
shown to avoid messing up the plot. The variational parame- 
ter used for the VCAn result was x = {V~\. The VCAn data 
was obtained for L = 6. For CPT as well as VCAn , we used 
a numerical broadening of 0"^ — 10~®. For comparison the 
results obtained by NRG (yellow) and FRG (dark brown) are 
shown^. 



semi-analytical calculation of Tk within VCAq . This 
additional factor is the same for all initial parameters 
(within the Kondo regime), it is particularly indepen- 
dent of A. If plotted over a scaled U axis, U' = iJ7, the 
VCAq result would lie on top of the NRG data. However 
using larger sizes of the interacting part of the reference 
system does not lead to much better results regarding 7. 
It is to be expected, that a significant improvement can 
only be obtained using exponentially large L. The CPT 
result shows a very different convergence behavior in L 
which is rather slow. 

An attempt was made to extrapolate the CPT data to 
i — > 00 by a simple scaling. It is interesting to ob- 
serve that this extrapolated curve coincides nicely with 
the VGA result (L = 6) in the low U region. However, 
we expect this extrapolation based on small L to be in- 
sufficient to capture the exponential scaling of the data 
in L. Note that since in CPT the self-energy is taken 
from the cluster, the CPT results for the effective mass 
coincide with ED results for systems of size L. 



2. Kondo spectral weight and half width 

Since the height of the Kondo resonance is fixed by 
the Friedel sum rule eq. the width and the weight 
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FIG. 14: (Color online) VCAn results for the spectral weight 
(blue) and full width at half maximum (olive) of the Kondo 
resonance as a function of interaction strength U . The vari- 
ational parameter used was x = {^}- A length of the in- 
teracting part of the reference system of L = 6 sites and a 
numerical broadening of 0^ — lO"*" were used for this calcu- 
lation. Data points marked with a circle were used for the fit 
of the exponential function in the Kondo region. The black 
line shows the Kondo temperature Tk as obtained by Bethe 
Ansatz calculations eq. (|24|) . 



(area) of the peak are proportional to the Kondo tem- 
perature Tk- Obtaining the spectral weight or FWHM 
of the Kondo resonance from the spectrum introduces a 
large uncertainty. Nevertheless we made an attempt, to 
get an idea of the behavior of Tk- We fixed the spectral 
weight by the first minimum to the left and to the right 
of the central peak (see also ref. Is^ ). In general the ef- 
fective mass and static spin susceptibility will yield more 
reliable results but it is instructive to compare these four 
ways of determining Tk- 

Shown in fig.[T3]is the evolution of the spectral weight and 
the FWHM of the Kondo resonance with increasing inter- 
action strength U . The data was acquired using VCAq 
with a variational parameter x = {V} for the particle- 
hole symmetric SIAM. Within the uncertainty, the same 
exponential behavior for the Kondo temperature Tk is 
obtained as by calculating the effective mass in VCAn . 



3. Static spin susceptibility 



The static spin susceptibility Xm is given by the linear 
response to an applied magnetic field B in z direction 



XrniU) 



dB 



(27) 



5=0 



In the Kondo regime this quantity too becomes inversely 
proportional to the Kondo temperature. For the calcu- 
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FIG. 15: (Color online) The static spin susceptibility Xm 
eq. (|27|l is shown as a function of interaction strength U. The 
variational parameter used was x — {V}. The data was ob- 
tained for L = 6 sites and 0^ = 10"*^. For comparison the 
results obtained by NRG (blue) and FRG (green) are showr>^. 



lations in this section we introduce an additional spin 
dependent term in the impurity Hamiltonian (eq. 



n 



B 



magnetic 



(28) 



The static spin susceptibility Xm as obtained with VCAq 
is shown in fig. 1151 The variational parameter used was 
x = {y}. As a reference the results of NRG and FRG^ 
are shown. The behavior of the VGA result is good for 
small interaction strength U. The VGA result shown for 
L = 6 appears already converged while the GPT result 
would require much larger systems. 

We would like to highlight that VGAo reproduces an 
energy scale Tk- Results from direct calculation of Tk, 
calculation of the effective mass m* and the static spin 
susceptibility Xm yield the correct functional form but 
not the right exponent. 



F. Benchmarking CPT/VCA against continuous 
time Quantum Monte Carlo 

In this section we compare GPT/ VGA results to QMG 
data. We obtained the Monte Garlo results using the 
continuous time Quantum Monte Garlo (GT-QMG) code 
of the TRIQS^ toolkit and its implementation of the hy- 
bridization expansion (GT-HYB)^ algorithm using Leg- 
endre polynomials^^. This method enables access to very 
low temperatures and is especially suited to obtain low 
energy properties^. The GT-QMG data provides statis- 
tically exact and reliable results to test our data. 
All GT-QMG calculations were done for a single impurity 
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orbital sdU = 0.8 and e/ = —0.4. We used a semicircular 
hybridization function with half bandwidth D = 2 and 
V = 0.3162. This setup corresponds to the same model 
under investigation here. The value for the interaction 
strength U = 0.8 was chosen because of the relatively low 
expected Kondo temperature of Pk = « 100. For 
all calculations 1.2 x 10^ MC updates where conducted, 
with a sweep size of 100 updates, plus a 10% thermaliza- 
tion period. 

To ensure that the Kondo resonance is correctly repro- 
duced by CT-QMC we evaluated the Matsubara Green's 
function for various values of inverse temperature (3. 
The height of the Kondo resonance is given by the 
Friedel sum rule eq. (f23| to be 9m(G//(icj„ = 0)) = — 10 
for the parameters used here (A = 0.1). To obtain 
3m(G'//(ia;„ = 0)) we extrapolate twice, first in iujn — )■ 
for each (3, then we use these results and extrapolate to 
T — >■ 0. The extrapolation to icj„ — > is done linearly us- 
ing the first two Matsubara frequencies. The imaginary 
part of Gff{iu!n) and the extrapolated value to iuj„ — )■ 
are shown in the inset of fig. [12] for /3 € [10, 1200]. Those 
extrapolated values are plotted as a function of tempera- 
ture (fig. [T5|) . These data points are then extrapolated to 
T — > using a fit by a rational model function. The re- 
sult clearly shows the onset of the Kondo resonance when 
the temperature is lowered below the Kondo temperature 
Tk- The extrapolation to T = shows very good agree- 
ment {^m{G f f{iuj„ = 0)) w —10.1) with the result ex- 
pected from the Friedel sum rule within the uncertainty. 
It is important to note that the CT-QMC results con- 
verge very nicely in /3. Although for higher /? lower Mat- 
subara frequencies become available, the overall shape of 
the Green's function does not change significantly. 
Therefore we may compare the T = CPT/VCA results 
for the Green's function and self-energy to the CT-QMC 
data. The Matsubara Green's functions of the impu- 
rity f-orbital Gffiiujn) obtained by CT-QMC (/3 = 400), 
CPT and VGA arc shown in fig.[I71 We use /3 = 400 
as a compromise between low temperatures and still re- 
liable CT-QMC results (within manageable computation 
time). The /? = 400 result was obtained using 65 Leg- 
endre coefficients. A detailed analysis has shown that 
this number is sufficient to get high frequency moments 
of the self-energy E accurately. The VCAn results were 
obtained with one variational parameter x = {V} for 
U = 0.8, A = 0.1 and 0+ = 10"*^ in the particle- 
hole symmetric case. For the CPT calculation we used 
the same parameters. For both methods we considered 
lengths of the interacting part of the reference system 
of L = 2,4,6,8 and 10. The VGA resuh lies near the 
CT-QMC data but underestimates the slope of the curve 
at low iojn- The VGA result provides a huge improve- 
ment upon CPT for the lengths of the interacting part 
of the reference system shown here. The real part of 
Gff{ioJn) is exactly zero within CPT/VCA as it is sup- 
posed to be. Note that the value of Gff{iujn = 0) which 
is fixed by the Friedel sum rule is exactly reproduced 
within CPT and VGA for the particle-hole symmetric 
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FIG. 16: (Color online) CT-QMC result for the imaginary 
part of the impurity Green's function extrapolated to iun ~ 0. 
An extrapolation to zero temperature is attempted, which 
yields a good agreement with the result predicted by the 
Friedel sum rule (green circle) within the uncertainty (red 
triangle). The inset shows the imaginary part of the impurity 
Green's function for various /? (see legend) and the extrapo- 
lated points at iujn = 0. 



case. The same is shown for the self-energy of the im- 
purity f-orbital I]//(ia;„) in fig. 1191 From the imaginary 
part of Tiff{iu)n) one can infer the convergence of the 
CPT/VCA result with larger length of the interacting 
part of the reference system L. The real part of the self- 
energy (3ffe(S//(ia;„) = fj, = —ef — ^ — 0.4) is again 
exactly reproduced within CPT/VCA. 
In the following, we discuss the self-energy Y,{iujn) for the 
two interesting cases of very low and very high Matsub- 
ara frequency. We start out by conducting an expansion 
of the self-energy E(z) for high Matsubara frequencies 
(z = iujn — > oo) which shall be outlined here briefly. The 
self-energy matrix is defined by 



E(z) 



G— 1 ~i 


z-T-G"- 



Here T is the one-particle part of the Hamiltonian. In 
the particle-hole symmetric case considered here it con- 
tains all the hoppings as well as the on-site energy of the 
impurity ef — — -j- We conduct a series expansion in 
powers of of E(z). Apart from the real constant Tu 
all z-dependent terms of 5]jj(z) are anti-symmetric in z. 
Therefore even powers in z^^' , / > vanish. Expanding 
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FIG. 17: (Color online) Comparison of the Matsubara impu- 
rity Green's function Gfjiiuin) obtained by CT-QMC (/3 = 
400), CPT and VCAn. The CPT/ VGA results were obtained 
for L = 2, 4, 6, 8 and 10. The real part shown in the lower 
part of the figure is zero. The Friedel sum rule prediction of 
Q'm(G//(ia;„ = 0)) = —10 is fulfilled by all methods. The leg- 
end of this figure serves as well as legend for fig.[T8]and fig. 1191 
That is why the last entry (large iujn exp. (see eq. (|29|l ')'l is 
displayed in the legend but is missing in the graph of this 
figure. 



the Green's function G(z) yields for the self-energy Ti{z) 

oo 

E(z) = -t-z^(-i)™a:'" , 

m— 1 

oo 

X = z " ■ C„ , 



(-1)" (*„|a](AH)"aJvI,o) 



where /S.'H = % — ljq and wq is the ground-state energy 
of i-L. Collecting powers of z yields a cumulant-like ex- 
pansion for the self-energy S(z) 

oo 

E(z) = ^z-"E„ , where 

n=l 

So = -T + Ci , and 
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FIG. 18: (Color online) Comparison of the self-energy of the 
impurity Tijjiiuin) times energy a;„ obtained by CT-QMC 
{13 = 400), CPT and VCAn. The CPT/ VGA resuhs were 
obtained for L = 2, 4, 6, 8 and 10. CPT as well as VCAn 
become exact for high Matsubara frequencies. An expansion 
of E(iaj„) for large iujn eq. (|29|l is additionally shown (straight 
line at — (y) )• The legend for this figure is the same as for 
fig. ll7l and is displayed there. 



where the self-energy at the impurity f-orbital is the 
only non- vanishing matrix element of E^j. This result 
is plotted as a reference in fig. [121 Due to the nature 
of the CPT/ VGA approximation these methods always 
yield the exact self-energy for high Matsubara frequency 
as shown in fig.fTSl 

The low energy properties examined in the previous 
section depend basically on the slope of the Matsubara 
Green's function at (iaj„) = 0+. The results shown in 
fig.lTTl and fig. [T^ show that this slope is underestimated 
by CPT/ VGA in comparison to CT-QMC, at least at 
the small lengths of the interacting part of the reference 
system available. 

The above results suggest a possible application of VGA 
as an impurity solver for zero temperature DMFT. 
The results would not suffer from a bath truncation 
error as in exact diagonalization based DMFT. A big 
advantage would be the low demand on computational 
power of VGA as well as the approximate reproduction 
of the main features of the local density of states (i.e. 
Kondo resonance and high energy incoherent part of the 
spectrum) . 



Here we consider the zeroth and first order in z ^ only 
and obtain for E(za;„) 



E//(iw„) = 




G. Introducing a symmetry breaking field 

We explore the possibility to improve the VGA results 
achieved by varying the internal single-particle param- 
eters of the model by introducing a symmetry breaking 
"spin flip field" at the impurity f-orbital. The term added 



18 




« 0.405 - 
s 

3= o-*- 

ir 0.395 - 

0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 
(0 

n 

FIG. 19: (Color online) Comparison of the imaginary part 
of the self-energy of the impurity ?sm(E//(ia;„)) obtained by 
CT-QMC (13 = 400), CPT and VCAn. The CPT/VCA re- 
sults were obtained for L = 2, 4, 6, 8 and 10. An expansion of 
'E{iu}ri) for large ia;„ eq. (|29[l is shown in addition (magenta 
line which diverges at zero). CPT/VCA always reproduces 
the exact self-energy for high Matsubara frequencies. The 
legend for this figure is the same as for fig.[T7]and is displayed 
there. 

to the impurity Hamiltonian cq. ([U 

Hmp = (/| 4 + /| /^) , (30) 

explicitly breaks the conservation of spin in the cluster 
solution. We are interested in the model with a physi- 
cal parameter = so this variable may only attain 
a finite value as a variational parameter B'^ in the refer- 
ence system. We investigate the particle-hole symmetric 
model a.t V = 0.3162 and t ~ 1. Our findings indicate 
that any finite value of B'^ splits the Kondo resonance 
and has thus to be discarded on physical grounds for the 
system under investigation. 

While this prevents the application of this field to im- 
prove the VGA results, it gives very nice insight in the 
physics of the SIAM as described by CPT/VCA. Wc find 
that a critical interaction strength Uc depending on the 
length of the interacting part of the reference system ex- 
ists which separates solutions which would prefer a finite 
B'^ from those which would prefer B'^ = 0. The critical 
interaction strength for i = 4 is given by Uc/A w 4.3. 
The grand potential fJ— Oq is plotted for various inter- 
action strengths U in fig-dOl For an analogous calculation 
for L = 6 site interacting clusters a value of Uc/A w 4.1 
is achieved. The mean field result would yield a critical 
interaction strength Uc/A = tt for the parameters used 
here. We interpret this value as a signature of the onset 
of local moment behavior. The values for Uc are of course 
not to be taken literally, they depend very much on the 
finite size of the cluster under investigation. 
The splitting of the Kondo resonance caused by a non- 
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FIG. 20: (Color online) Grand potential - (eq. (f20)) ') 

as a function of the interaction strength U/A (see legend). 
The data was obtained by studying a L = 4 site interact- 
ing cluster coupled to a semi-infinite lead. The numerical 
broadening used was 0"*" — 10~^. The crosses indicate the 
respective minimum of the grand potential. There exists a 
critical Uc/A. « 4.3 above which a finite B'^ is preferred by 
the system. 

zero variational field B'^ is shown in fig. [211 The value 
of J7/A = 12 used for this calculation lies in the region 
above Uc where the system prefers a nonzero field B'^. 



VI. CONCLUSIONS 

In this work we have applied the variational cluster 
approach (VCA) to the single impurity Anderson model. 
We devised a cluster tiling applicable to this non- 
translationally invariant model which leads to a cluster 
with a discrete spectrum and an environment having a 
continuous spectrum. Wc have derived an expression for 
the change of the grand potential originating from the 
coupling of the impurity to the semi-infinite bath. 
We have compared results for the single-particle dy- 
namics to data obtained by exact diagonalization and 
cluster perturbation theory (CPT). We found that the 
variational extension made by the VCA is vital for a 
good reproduction of the expected behavior of the SIAM. 
The CPT/VCA spectra both yield a Kondo resonance 
in the impurity density of states with the correct height 
as predicted by the Friedcl sum rule. A close look at 
the Kondo resonance shows that the VCA is able to 
reproduce the resonance and the functional behavior 
for the Kondo temperature in a remarkable way. The 
Kondo temperature is expected to show exponential 
behavior in interaction strength in the Kondo regime. 
VCA yielding an exponential behavior however tends to 
underestimate the exponent. Comparison of dynamic 
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FIG. 21: (Color online) The splitting of the Kondo resonance 
caused by an applied magnetic field in x direction is shown 
for different values of the auxiliary field B'^. The plots were 
obtained using VGA (i.e. the physical field Bx is always zero). 
Instead of taking the parameter B'^ at the stationary point of 
the grand potential (this value would be B'^/ A ~ 1.9 for the 
parameters used) we explicitly plug in a fixed value for B'^. 
The length of the interacting part of the reference system 
used was L — 6 for the model parameters U /A — 12. The 
numerical broadening used was 0^ = 10~^. 



quantities to continuous time Quantum Monte Carlo 
solidifies the origins of this behavior. The high energy 
incolierent part of the spectrum sliovifs strong finite size 
effects within CPT which are partly removed by virtue 
of the VGA. VGA furthermore reproduces the expected 
position and width of the high energy part of the 
spectrum. For the asymmetric model, the Friedel sum 
rule is fulfilled in all parameter regions implying that 
the Kondo resonance is pinned at the chemical potential 
in the Kondo region. In addition a self consistent 
formulation of the VGA, previously introduced in the 
context of non equilibrium problems^, was explored. 
Results obtained by the self consistent approach show 
agreement with results obtained by VGA based on the 
grand potential for the density of states of the impurity 
f-orbital. Thereby the positions of the spectral features 
agree very well with the traditional VGA result, while 
the spectral weight distribution may deviate especially 
for small values of interaction strength. Comparison 
to results obtained from Bethe Ansatz, renormalization 
group approaches and data obtained from X-Operator 
based calculations show reasonable agreement for all 
quantities investigated. 

In conclusion, while there are certainly more accurate 
methods to deal with a single quantum impurity model^^, 
especially at low energies, our work shows that VGA is a 
fiexible and versatile method which provides reasonably 
accurate results with modest computational resources. 
Here, the VGA self consistency condition proves to be 
crucial. This allows to obtain the same accuracy that 



GPT would provide with a much larger, inaccessible, 
interacting cluster size. One of the advantages is the 
ficxibility of the method, i.e. it is straightforward to ex- 
tend it to many impurities, non equilibrium problems^, 
etc. In the spirit of NRG^ one could improve on the 
present results by carrying out an appropriate unitary 
transformation on the bath such that the bath hoppings 
decay with increasing distance from the impurity. In 
this work, we do not include such an improvement and 
choose a constant hopping sequence (eq. since our 
goal is to benchmark VGA/GPT only. Nevertheless, a 
hybrid approach combining NRG and VGA would be an 
interesting extension of the present work. 
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Appendix A: Grand potential 



Here we outline the proof of eq. p9|) . We start out 
from eq. ([TS]), i.e. 



An = -Tr In U - Tg 



Taylor expansion yields 





'(Tg)" 


+ Tr 


'(Tg)" 
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n [ 




CC 







(Al) 



Due to Tee = each term ggg in the first trace occurs 
only in the form 

Sec Tceg'eeTec • 

The expressions in the second part of eq. (jAip can be 
modified by a cyclic permutation of the factors in the 
argument of the trace 

Tl- (Tecgcc • • • Tcegee) === Tr (gee ' ' ' TeegeeTec) , 

in which g^^ again only occurs in the dressed form gee. If 
we replace all occurrences of TeegeeTec by gee then there 
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are no matrices Tee and Tec respectively, left. Hence we 
can as well introduce in eq. ([TS]) the following replace- 
ments 

gee ice : Tee ^ lee , Tec ^ Ice , H ^ Icc , 

as it leads in the series expansion to the same expressions. 
The argument in eq. (|18p then assumes the form 



1 



with the abbreviation 5 := IL — Teegcc- Prompted by 
the Schur complement, the matrix can be factorized into 
upper and lower triangular block matrices 
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-^'-'ge 
1 



such that the determinant is easily computed, since the 
determinant of the second matrix is 1 and the determi- 
nant of the first matrix is simply the product of the de- 
terminants of the diagonal blocks, resulting in 



det 



b -i 



det (6) det (U 



det 6 - ge 



= det - (Tee + TeegeeTee)ge 

The final result for eq. reads 



An = -Tr In ( lee - Seegec 



Appendix B: Behavior of tlie Friedel sum rule 
within ED, CPT and VGA 

It is possible to gain a somewhat deeper understanding 
of the behavior of the FSR within ED/CPT/VCA by 
considering the local Green's function at the impurity 
f-orbital 



G//(z) = (z-e/-r(z)-S](z))- 



(Bl) 



where T(z) is the contribution due to the single-particle 
terms of the hybridization and the self-energy due 
to the local interaction (for details see app. [C|. In the 
following we consider the particle-hole symmetric case. 
We are interested in the behavior of the retarded Green's 
function G^°^(a;) = G//(a; -I- iO~^) at the Fermi energy 
(id = 0); which we investigate by taking the limit on 
the Matsubara axis Gf}{0)^ lim G//(w). We expand 

the self-energy E(iz/) up to linear order in v and rewrite 
the expression using the definition of the effective quasi 
particle mass m* (see eq. ([26])) 



E(w) ~ -£/ + j3m(E(0)) + iv 
Ri —ef + w(l — m*) . 



(99m(E(w)) 



0+ 



Inserting into eq. (|Bip we obtain 
Gf}iO) = ^hm Gffiz,,) 



= lim 



i [m*iy - 3m(r(i:/)] - ^c{T{iiy) 



(B2) 



We will now investigate two separate, general cases of an 
ED and a CPT/ VGA treatment of the Green's function. 
From eq. (|B2p it is easy to sec that the remnant to* of the 
impurity self-energy and therefore the self-energy itself 
does not contribute to the Friedel sum rule. In outlining 
how to notice this, we simultaneously show that a dis- 
crete spectrum of the conduction band (as obtained for 
example in ED) will not fulfill the Fridel sum rule. Con- 
sider an arbitrary discrete spectrum of the conduction 
electrons with hybridization 



^ ' 7 7/ — /, 



with excitation energies w^. Splitting into real and imag- 
inary parts and inserting into cq. (jB2p gives 

Qm{Gf}{0))= liin^ 

- (to* + V'^A{iy)) V 



(m*2 + 2m*V'^A{v) -t- V^A^vf) v"^ + V^B[vY 



with 



A* 9 I ? 



and 



Upon neglecting the weak dependence of Aiv) and B{y) 
on V one obtains, in this case, in the limit — ^ 0: 
3m(G7;(0)) ^ if all cj^ 7^ or 5m(G7/(0)) -00 if 
any = 0. We would like to further illustrate this for 
the specific model considered in this work (i.e. an impu- 
rity coupled to a semi-infinite chain with open boundary 
conditions) 



L-l 



i=0 



where we supressed spin indices. Wc obtain by using 
the equation of motion for the local impurity f-orbital 
Green's function 



>> 



= ([c/,4]+)- << [c^.,-H]_;4 >> , 
for an even number of sites (including the impurity) 

3m(rcvon(«i^)) = RL{v,t) , 

V 
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and for an odd number of sites 

9m(roddM) = -V^iyRL{'^,t) , 

where Rl{v, t) is a rational function which is weU behaved 
upon taking the hmit — > 0+ (i.e. it approaches a con- 
stant Um Ri^iv^t) = fLit)) a-iid the real part is always 

zero. Upon inserting into eq. (jB2p one can easily verify 
that for an even number of sites this yields zero spectral 
weight at the Fermi energy, while for an odd number of 
sites it yields —oo (i.e. there is a pole exactly at a; = 0). 
Therefore the ED results alternate with even/odd system 
size between pf™(0) = and pf'^iO) = oo. This result 
shows that the FSR is always violated in ED because a 
finite value of the impurity density of states at the Fermi 
energy may only be obtained using an artificial numerical 
broadening. It furthermore shows that all terms involv- 
ing m* go to zero and cannot contribute to the sum rule. 
Now we turn to the case of CPT/VCA, where the conduc- 
tion electron hybridization takes the form (see eq. ([TT 



I , 



for the model considered in eq. ([T]) (i.e. a semi-circular 
density of states of the conduction electrons). The 
CPT/VCA Green's function of the physical system is 
then given upon insertion of this F into eq. ()B2p 

5m(G7.*(0)) = lim 



2t^ 



Which yields upon expansion of the square root up to 
linear order in v 



T 



and when the limit — > is taken 
3m(G7;(0)) 

which is exactly the value predicted by the FSR (eq. ([231)) 



3m(G7;(0)) 



1 . TT 2 

-TT — - sin ( — ) 



1 

"a 



This result is independent of the size of the reference 
system. Away from particle-hole symmetriy (occupation 
of the f-orbital not one) the calculation becomes more 
tedious. The numerical VCA calculations however show 
that a pinning of the Kondo resonance at the Fermi 
energy is obtained in contrast to CPT at small L. 



Appendix C: Semi-analytical expressions for VCA of 
the two-site problem 

To gain a better understanding of the behavior of 
the low-energy properties of the SIAM we solve a small 
system semi-analytically. We obtain the scaling of the 
Kondo temperature Tk with interaction strength U 
within VCAn as well as VCAsc • A reference system 
consisting of a two site cluster and an infinite environ- 
ment is used. 

The Kondo scale will be determined from the effective 
mass eq. (|26p using the self energy of a two site cluster 
eq. (jC3p which leads to 



m*{U) = 1 + 



1 

36 



(CI) 



The Kondo scale Tx is inversely proportional to m*{U) 
eq. psp . Therefore within this approximation the behav- 
ior of the optimal cluster parameter V' governs the low 
energy physics. 

To determine the optimal hopping V' the Green's func- 
tion of the reference system g is calculated 

/2-E'(z) V 
g-\z) ^ \ V z 

V Q'-\z)^ 

The CPT/VCA Green's function G 

lz-Yl{z) V 
G-i(z) V z t 

\ t g;;1(z), 

is obtained by eq. ([6]) using 

Ay 
T = I Ay 

-t 

Here Gn = G// eq. (|C2p corresponds to the impurity f- 
orbital, G22 = Gss the second site in the interacting part 
of the reference system and G33 = Gee eq. (fTO|) the semi- 
infinite environment. The cluster parameter V' is given 
by the physical parameter V plus the variation Ay. Note 
that in here we work with the reduced expressions for f2 
and G justified in app. |^ Schonhammer and Brenig cal- 
culated the Green's function of the correlated orbital for 
this model perturbativeley and showed that their expres- 
sion becomes exact in the limit of vanishing bandwidth^!. 
This is exactly the case considered here, where the im- 
purity f-orbital is coupled to a single non interacting site 
providing a bath with vanishing bandwidth. They ob- 
tained 



1 



2-r'(z)-i]'(z) 

where the hybridization r'(z) in our case is given by 

r'(^) = — , 

z 



(C2) 
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and the self-energy "S'iz) is given by 



Evaluation of this expectation value in the cluster yields 



S'(z) = 



4 



z - 9r'{z) 



(C3) 



From this all elements of the cluster Green's function may 
be obtained by the equation of motion technique. 
To be able to calculate the grand potential il, the ground 
state energy of the interacting part of the reference sys- 
tem, needs to be obtained (which can be done for 
example by diagonalization of the Hamiltonian matrix 
or by an integral over the Green's function) 



Using the expression for the CPT/VCA Green's function 
g and the ground state energy and taking the derivative 
of the grand potential cq. ([20]) with respect to AV one 
is able to obtain an integral expression which allows to 
determine V within VCAn : 



dn{AV) ^ 1 



l^T 1^ (a - T(Ay)g(iw, Av)r^ \^ 
(VAyT(Ay)g(zt^,Ay)) 
+ (T(Ay)VAyg(iu;, AV)) ) ) ) - . (C4) 



The resulting V'{U) is shown in fig.[21](left) and is used 
to calculate the effective mass eq. (jCip shown right in 
the figure. The effective mass shows exponential behav- 
ior but the exponent does not match the Bcthe Ansatz 
result as discussed in sec. IV El 

Next we attempt to obtain the VCAsc solution for the 
two-site problem. The only variational parameter is 
V and therefore we determine the expectation value of 
J2 < ft^ia > self consistently. Here 1 denotes the 

impurity's s-orbital. Since we are considering a spin- 
symmetric model we sum over both spin directions and 
denote this expectation value as < /^c > in the following. 
The hopping expectation value is given by 



n 



did G/c(ii^) 



2 f°° 

< Pc >clustcr= / duj 

n 



4{AV + V) {9{AV + V)^ + w^) 
36{AV + V)^ + (C/2 + 40(Ay + V)^) + Aw^ 



0.35 
0.3 
0.25 



semi-analytic VCA^, L=2 

— semi-analytic VCA^^, L^2 
+ VGA U2 (from numerics) io* 




semi-analytic VCA^, L=2 

semi-analytic VCA^^, L=2 
o NRG 
Bethe Ansatz 











FIG. 22: (Color online) Left; Optimal parameter V' of the 
reference system as obtained by the semi-analytical equations 
for VCAn eq. (fC4l) and VCAsc eq. (fCS)) . As a reference the 
L = 2 data of our numerical simulation is shown too. Right: 
The effective mass eq. (|C1|) obtained by the optimized param- 
eter V' of the reference system (see left figure). Additionally 
shown is the Bethe-Ansata^ eq. (|24|) and NRG^ result as a 
reference. 



Evaluation of this expectation value in the total system 
gives 



< .fc >CPT 



du 



_j_ w(U^+i6(AV+VY+iw'^)[w+y^Ii^T^) 



9(AV-l-y)2+u 



Upon requiring the two expectation values to coincide 



< fc >clustcr = < Pc >CPT , 



(C5) 



the optimal value of AV is obtained numerically. The 
resulting V'{U) is plotted in fig. [22] (left) and is used to 
calculate the effective mass eq. (jCl|) shown right in the 
figure. The effective mass does not show exponential be- 
havior in VCAsc ■ 
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